Crack propagation as a free boundary problem 
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A newly developed sharp interface model describes crack propagation by a phase transition pro- 
cess. We solve this free boundary problem numerically and obtain steady state solutions with a 
self-consistently selected propagation velocity and shape of the crack, provided that elastodynamic 
effects are taken into account. Also, we find a saturation of the steady state crack velocity below 
the Rayleigh speed, tip blunting with increasing driving force and a tip splitting instability above a 
critical driving force. 

PACS numbers: 62.20.Mk, 46.15.-x, 46.50.+a, 47.54.-r 
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One of the most challenging riddles in physics and tech- 
nology is the phenomenon of fracture, as it gives rise 
to material failure on all scales. Most fundamentally, 
the initiation of crack growth is due to a competition of 
the release of clastic energy and an increase of surface 
energy, which has been pointed out by Griffith and 
has been used to describe many features of cracks 0- 
The interpretation of brittle fracture as the successive 
breaking of atomic bonds is in agreement with models of 
sharp crack tips, but still the theoretical predictions de- 
pend significantly on empirical interaction potentials Q . 
On the other hand, if dissipative plastic effects or large 
scale deformations are important, the crack tips are ex- 
tended and have a finite tip radius r$. Here, a detailed 
description of fracture necessitates equations of motion 
for each interface point of the extended crack instead of 
just the mentioned integral energy balance. A full mod- 
eling of fracture should then not only predict the growth 
velocity but also determine the entire shape of the crack 
sclf-consistcntly. It was proposed that the characteristic 
length scale of the crack tip tq is selected by the threshold 
of plastic yielding pj. However, approaches of this type 
require the introduction of dynamic theories of plasticity 
which are usually much more speculative and less ver- 
ified than the ordinary linear theory of elasticity. The 
lack of suitable equations of motion becomes apparent 
also in the regime of fast fracture, where the experimen- 
tally observed maximum crack speeds are far lower than 
the theoretically expected Rayleigh speed 0, • Beyond 
a critical velocity, an unpredictable tip splitting of the 
crack can occur, producing oscillations of the crack speed 



Specific equations of motion have been implemented in 
various phase field descriptions 0, B IE ilil III! ■ They 
provide a stimulating approach to describe fracture as a 
moving boundary problem and go beyond discrete mod- 
els. However, additional parameters are introduced in 
these models in comparison to a conventional linear elas- 
tic theory, and the scale of the patterns in the tip region is 
therefore typically selected by numerical parameters like 
the phase field interface width. A purely static elastic 
description together with macroscopic equations of mo- 



tion [l^ provides a well defined sharp interface limit, but 
suffers from inherent finite time singularities in this limit, 
which do not allow steady state crack growth. Based on 
these insights we recently developed a continuum theory 
of fracture which resolves t his p roblem by the inclusion 
of elastodynamic effects 0, Il4j . The main advantage of 
this model is that it relies only on well established ther- 
modynamical concepts. Since in the phase field descrip- 
tion |bj] an extended hierarchy of length scales has to be 
resolved, expensive large-scale simulations are necessary 
to predict the steady state growth of cracks: The system 
size must be much larger than the crack tip scale, which 
itself must exceed the phase field interface width and 
the numerical lattice parameter significantly. Further- 
more, only after a long relaxation time this fully time- 
dependent description can converge towards a steady 
state solution. 

Here we propose an alternative approach which is 
specifically dedicated to the fast steady state growth of 
cracks. The limit of fully separated length scales is per- 
formed analytically, leading to a very efficient numerical 
scheme. Furthermore, this approach reaches the valid 
steady state regime with a self-consistently selected crack 
shape much faster, which also accelerates the computa- 
tions dramatically. The method is based on a multipole 
expansion technique, and the satisfaction of the elastic 
boundary conditions on the crack contour reduces to a 
linear matrix problem, whereas the bulk equations of dy- 
namical elasticity are automatically satisfied. Neverthe- 
less, finding the correct crack shape and speed remains a 
difficult nonlinear and nonlocal problem. In this Letter, 
we formulate the free boundary problem of crack propa- 
gation based on a phase transformation model, and solve 
it numerically. 

Continuum Model of Fracture. Imagine that the crack 
is filled with a soft condensed phase instead of vacuum, 
and the growth is then interpreted as a first order phase 
transformation of the hard solid matrix to this soft phase 
0i Q ■ The inner phase becomes stress free if its 
elastic constants vanish. For simplicity, we assume the 
mass density p of both phases to be equal; together with 
the premise of coherency at the interface this implies that 
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the solid matrix is free of normal and shear stresses at 
the crack contour, i.e. o~ nn = a nT = 0, which serves as 
boundary conditions. In the bulk, the elastic displace- 
ments Ui have to fulfill Newton's equation of motion, 



= pui. 



(1) 



The difference in the chemical potentials between the two 
phases at the interface is given by 



A^i = ( ->i-rr,,(,, - 7K), 



(2) 



with 7 being the interfacial energy per unit area; the in- 
terface curvature k is positive if the crack shape is convex; 
the atomic volume Q appears since the chemical poten- 
tial is denned as free energy per particle. For elastically 
induced phase transitions the motion of the interface is 
locally expressed by the normal velocity 
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(3) 



with a kinetic coefficient D with dimension [D]=m 2 s _1 . 
We have investigated this minimum model of fracture 
given by Eqs. Q - © by phase field modeling in Q and 
have demonstrated that it describes fast crack propaga- 
tion provided that dynamical elasticity is taken into ac- 
count. 

The Multipole Expansion Method. We discuss the 
steady state propagation of a semi-infinite crack in an 
isotropic medium. We assume a two-dimensional plane 
strain situation and mode I loading, which means that 
the applied tensile forces act perpendicularly to the crack. 
The negative x axis is aligned along the semi-infinite 
crack, see Fig. We start with the description in the 
laboratory frame of reference and assume that the crack 
is moving with a given constant velocity v. Following 
Ref. we introduce two real functions </>(x, y, t) and 

ip{x, y, t) which are related to the displacements Uj as fol- 
lows, 

d<f> dip d<p dip 

x dx dy ' v dy dx 

With this decomposition, the bulk equation JIJ decouples 
to two wave equations, 

c 2 V 2 = Ofa, c 2 s V 2 iP = d^. (4) 



Here, c d = yjE{\ - v)/p{\ - 2v){\ + v) and c s = 
y/E/2p{l + v) are the dilatational and shear sound speed 
respectively, E is Young's modulus and v the Poisson ra- 
tio. 

In a steady state situation the time derivatives in 
Eqs. Q vanish in a co- moving frame of reference (x — > 
x — vt) and they become Laplace equations there, 

d 2 cf> d 2 (p _ d 2 ip d 2 ip 
dx 2 dy 2 , ' dx 2 dy 2 
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(5) 



We have introduced rescaled coordinates perpendicular 
to the crack, y d = a d y and y s = a s y. with a 2 = l — v 2 /c 2 
and a 2 = 1 — v 2 /c 2 . For a straight crack with a sharp 
tip, the solution obeying mode I symmetry and the usual 
a ~ r -1 / 2 behavior is 



a 3/2 3 

A r d < cos-i 



ili = -B rl/ 2 sin -6 



in rescaled polar coordinates which are related to the 
co-moving cartesian coordinates via x — r d cos 9 d = 
r s cos 9 S , y d = r d sin 6 d and y s = r s sin 9 S . For this 
mode, the boundary conditions on the straight cut and 
the matching to the far field behavior demand 



A = 



8(1 +!/)(! + a 2 s ) 
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K dyn , (6) 
(7) 



where K dyn is the dynamical mode I stress intensity fac- 
tor 0. 

In order to solve the elastodynamic problem of a crack 
with finite tip radius ro, we use a multipole expansion, 
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Aq cos 3#d/2 + T cos U" 

n—l a 



„ r 3/2 



B sin 36 s /2 
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Each eigenmode satisfies the elastodynamic equations 
(JSJ. On macroscopic distances r from the tip, i.e. ro <C r, 
the crack still looks like the semi-infinite mathemati- 
cal cut and therefore exhibits the same far-field behav- 
ior. Thus, the coefficients an d Bo are determined by 
Eqs. 10 and Q, whereas all other modes decay fast and 
do not contribute to the asymptotics. Consequently, we 
obtain the formal stress field expansion, 



K, 



JV=oo 



dyn 



(2 7 rr) 1 /2 



(n)' 



where f^\{9 dl v) and f\^\{9 s ,v) are the universal an- 
gular distributions for the dilatational and shear contri- 
butions which also depend on the propagation velocity. 
The unknown coefficients of expansion can be found by 
solving the linear problem of fulfilling the boundary con- 
ditions a nn = o~ nr = on the crack contour. The tan- 
gential stress o~ TT is determined only through the solution 
of the elastic problem, and enters into the equations of 
motion and Q , leading to a complicated coupled and 
nonlocal problem. 

The strategy of solution of the problem is as follows: 
first, for a given guessed initial crack shape and veloc- 
ity, we determine the unknown coefficients A n and B n 
from the boundary conditions. Second, we calculate the 
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chemical potential and the normal velocity at each point 
of the interface. Afterwards, the new shape is obtained 
by advancing the crack according to the local interface 
velocities. This procedure is repeated until the shape of 
the crack in the co-moving frame of reference remains 
unchanged. It provides a natural way to solve the prob- 
lem, as it follows the physical configurations to reach the 
steady state. Originally, this idea was developed in the 
context of dendritic growth [J^ • Then the following re- 
lation between the local normal velocity and the steady 
state tip velocity v holds: 

v n — v cos r\ = 0, (8) 

with rj being the angle between the normal to the crack 
contour and the x axis. Alternatively to this "quasi- 
dynamical approach" , we also directly solve the nonlinear 
equation (jHJ as a functional of the crack shape and the tip 
velocity v by Newton's method complemented by Pow- 
ell's hybrid method [Tsl flf| ("steady state approach"). 

Our results are obtained with a finite number of modes 
N, and we performed the extrapolation N — > oo and 
found only minor deviations of a few percent from the 
results for N = 12 presented here. An extended discus- 
sion of the numerical details will be published elsewhere. 

Results and Discussion. We discuss crack growth in 
a strip geometry, with the width of the strip being very 
large in comparison to the crack tip scale. We introduce 
the dimensionless driving force A = K 2 tat (l — v 2 )/2E~f, 
where the static stress intensity factor K sta t is assumed to 
be given and finite. The relation between the static and 
the dynamic stress intensity factor can be easily found 
from energy considerations for the strip geometry in the 
spirit of |20j 

K dyn = Kstat (1 - v) — - — 5V -^- . (9) 

The multipole expansion technique confirms that the 
simple phase transformation model, based on Eqs. 

provides a selection of the steady state shape and 
the velocityof the crack. A typical crack shape is shown 
in Fig. n |2lJ. The dimensionless crack velocity v/vr 
(vr is the Rayleigh speed) and the dimensionless crack 
opening v^h/D as function of the driving force A are 
shown in Figs. [2] and [3] The simulations have been per- 
formed with Poisson ratio v = 1/3, which is the only 
remaining parameter. All results are obtained both by 
the "steady state approach" and the "quasi-dynamical 
code", and they are in excellent agreement with each 
other. Above the Griffith point A > 1 the shape of the 
crack obeys the equation — vy' = Dy" in the tail region 
x — > — oo, because the elastic stresses have decayed there 
[lij . Its general solution 

y(x) = t;+ b cxp(~vx/D) (10) 
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FIG. 1: Shape of the crack obtained for A = 1.3; h is the tail 
opening. 
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FIG. 2: Steady state velocity as a function of the dimension- 
less driving force A. The solid line corresponds to the steady 
state code, the squares to the quasi-dynamical code. Below 
the point A c ~ 1.14 the dissipation-free solution is selected 
by a microscopic length scale. Also for A < 1 the tip scale 
is not selected, and the presented solution is obtained for the 
specific parameter vrIi/D — 10. 



contains a growing exponential which is excluded by the 
boundary condition of straightness. y' = 0, and finally 
leads to the selection of both the steady state propaga- 
tion velocity and the crack openin g h (For more detailed 
counting arguments, see Refs. Il^jh 

In principle, one would expect steady state solutions 
for crack growth to exist for all driving forces beyond the 
Griffith point, A > 1. However, in the framework of the 
model, they exist only in the interval A = 1.14 — 2.5. 
At the limiting value, A ks 2.5, the propagation velocity 
tends to zero and the length scale h diverges. Neverthe- 
less, at this point the product vh/D remains finite, as 
it is required for finite driving forces. This termination 
of the steady state solution is surprising, as one would 
expect the tip blunting to continue to arbitrarily large 
values of A. 

At the lower limit, A r* 1.14. the steady state crack 
velocity is finite, but the tail opening tends to zero in 
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the framework of the model. At this point, the dissipa- 
tion becomes zero, and all energy is converted to surface 
energy apart from kinetic contributions which are trans- 
ported through the soft phase and out of the system. 
Below this point, i.e. for 1 < A < 1.14, dissipative solu- 
tions do not exist. Naturally, the tip scale should then 
be determined by an intrinsic microscopic length scale 
which is not contained in the present model. If it was in- 
troduced here explicitly, the behavior of the crack speed 
would behave as depicted by the dotted curve in Fig. 
then it would become zero at the Griffith point A = 1. 
Precisely this behavior near the Griffith threshold was 
observed in phase field simulations 0] , where this cutoff 
naturally appears as the phase field interface width. 

As we have already noted, the length scale of the crack 
tip becomes large for high driving forces, and therefore 
at least in this region our macroscopic theory should be 
valid. On the other hand, the velocity decays in this 
regime with increasing driving force, which is an counter- 
intuitive outcome of the model; nevertheless, the product 
vh/D which controls dissipation is monotonically grow- 
ing. Physically, it means that the dissipation is mainly 
increased due to tip blunting instead of a rise of the crack 
speed. 

We suppose that the solutions become unstable against 
a secondary Asaro-Tillcr-Grinfeld instability [2*3 beyond 
the point A ps 1.8, in agreement with previous conjec- 
tures ^3 and phase field simulations ^4|- This is indi- 
cated here by the change of sign of the tip curvature kq 
as shown in Fig. QJ corresponding to a tip splitting struc- 
ture. Since by construction we confine our investigations 
to symmetrical steady state solutions, we cannot capture 
the full asymmetric scenario. If we released this con- 
straint, the competition between the emerging new tips 
would lead to complicated non-stationary growth scenar- 
ios, as we have seen in phase field simulations 0] . 

Finally, we discuss the healing of cracks below the Grif- 
fith point, A < 1; the velocity v of the crack becomes 
negative in this regime. In contrast to the case of growth, 
one expects these steady state solutions to exist for arbi- 
trarily prescribed openings h and only the velocity v to 
be selected. This corresponds to the fact that the grow- 
ing exponential in Eq. (|l()(l automatically vanishes in the 
tail region. This prediction is numerically confirmed by 
our simulations, see Fig. [5] We note that without elastic 
stresses, i.e. for A = 0, the problem has a simple analyt- 
ical solution: x{y)/h = (l/7r) lncos(7ry/ft.) with velocity 
v = nD/h. 

Im summary, we have presented a continuum theory of 
fracture based only on the linear theory of elasticity and 
a phase transformation process at the crack surface; we 
employ a sharp interface method to find steady state solu- 
tions of crack growth and are able to predict the growth 
velocity and the self-consistently selected shape of the 
crack. Beyond a critical driving force a negative tip cur- 
vature indicates the transition to a tip splitting regime. 




FIG. 3: The dimensionless tail opening vrH/D as a function 
of the dimensionless driving force A. The solid curve cor- 
responds to the steady state code, squares belong to quasi- 
dynamical calculations. 
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FIG. 4: Tip curvature no as a function of the driving force A. 
The solid curve corresponds to the steady state code, squares 
belong to quasi-dynamical calculations. 



The results are in qualitative agreement with phase field 
simulations flij ]. 
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